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ABSTRACT 


The autoregressive (AR) model of a random process is interpreted in the light of the Prony’s relation 
which relates a complex conjugate pair of poles of the AR process in the z-plane (or the z domain) 
on the one hand, to the complex frequency of one complex harmonic function in the time domain 
on the other. Thus the AR model of a time series is one that models the time series as a linear com- 
bination of complex harmonic functions, which include pure sinusoids and real exponentials as 
special cases. An AR model is completely determined by its z-domain pole configuration. The 
maximum-entropy/autoregressive (ME/AR) spectrum, defined on the unit circle of the z-plane (or 
the frequency domain), is nothing but a convenient, but ambiguous visual representation. It is 
asserted that the position and shape of a spectral peak is determined by the corresponding complex 
frequency, and the height of the spectral peak contains little information about the complex ampli- 
tude of the complex harmonic functions. 
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INTRODUCTION 


Maximum-entropy (ME) spectral analysis, as the name suggests, seeks to maximize the infor- 
mation entropy of a random process. The method was originally proposed by Burg (1967). Owing 
to its superior resolution to the conventional Fourier spectral analysis, ME spectrum has since been 
extensively used in research fields such as speech processing, astronomy, and geophysics. Van den 
Bos (1971) later showed that the ME spectrum is formally identical to the autoregressive (AR) spec- 
trum, and hence that the ME spectral analysis is equivalent to the AR modeling of a random process 
(see e.g. Makhoul, 1975). Since the latter is much more succinct in concept, we shall in this study 
consider the ME spectral analysis from the AR point of view, and speak of the ME/AR spectrum. 

Suppose we have an equally spaced, real-valued time series |x(n); n=l,2,. . ,,N J, which, as usual, 
is a realization of a linear random process. The AR model of the time series is the following: 

M 

x(n)= 2 S k x (n-k) + a(n), n = M+l,. . ,,N (1) 

k=l 

where M is the order of the AR model, |S k ;k=l, 2,. . .,M} are called the AR coefficients, and 
Ja(n); n=l,2,. . .,N} is a realization of a zero-mean, white random process with variance <a 2 >. In 
this paper we are mainly concerned with the ME/AR modeling of harmonic processes. An important 
special case is given by a time series generated by an impulsive excitation a(0) at time zero while all 
subsequent value of a(n) are zero. Then x(n) (normally corrupted by some additive noise) gives the 
impulse response of the physical system, and equation (1) represents a truncated complex harmonic 
process (as we shall see later). In any event, the parameters M and { S k } completely determine the 
AR process. In practice, they are to be estimated from |x(n); x=l,2,. . .,N } by minimizing <a 2 >. 
Equation (1) can then be used (if so wish) in the prediction of future values of the time series by 
setting a(n)=0. In fact, (1) represents one important class of models in the linear prediction theory 
(see e.g. Box & Jenkins, 1970, Makhoul, 1975), and <a 2 > is often called the prediction-error 
power. 
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The ME/AR spectrum can be obtained from (1) by taking z transforms. Thus, 

X(z) = z M A(z)/H(z) (2) 

where A(z) is the z transform of a(n), and 

H(z) = Z M -Sj z M "l-. . ,-S M (3) 

is a polynominal of degree M in complex variable z. The spectrum, therefore, is 

P(z) = | Z 2M| <a 2 >/|H(z)| 2 (4) 

evaluated on the unit circle of the complex z-plane, CQ:z=exp(ico); or, in a more familiar form, 

M 

P(w) =<a 2 >/|l - 2 Sfc exp (-icok)l 2 (5) 

k =1 

for 0<co<27 r, where oj is the argument of z, or the angular frequency. 

There has been a number of methods proposed for the estimation of j j and <a 2 > (see e.g. 
Makhoul, 1975, Ulrych & Bishop, 1975, Ulrych & Clayton, 1976). They are all based on some sort 
of linear least-squares procedures. The determination of M, on the other hand, is a problem much 
more complicated and controversial. However, we are not concerned with any of the estimation 
methods in this study. In fact, our present concern starts only after the model parameters have been 
determined or given. Essentially we shall address the following question: What is the meaning of the 
AR model (1) and the ME/AR spectrum (5)? 

The AR Model and Prony’s Relation 

First let us concentrate on the polynominal H(z). According to equation (2), the zeros of H(z) 
arc also the poles of the AR process (1). In fact, AR models are often called all-pole models in engi- 
neering literatures. Since the polynominal coefficients 1,-S^ ,. . .,-Sjq are ad real (see equation 1), 
the fundamental theorem of algebra ensures that these poles are either real or in complex conjugate 
pairs. Suppose, for the time being, that M is even and all the poles come in complex conjugate pairs. 
Let M=2Mj . Then in general we can express H(z) as 
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( 6 ) 


M t 

H(z) = II [z-exp(iaj)] [z - exp(-ioj*)] 
j=l 

with poles at {exp(iaj), exp(iOj*); j= 1 ,2,. . .,Mj } . The meaning of the complex quantities { aj ; j=l ,2, 

. . .,Mj j becomes obvious once we see that the AR time series (1), with a(n)=0, can be expressed as 
the following linear combination of Mj complex harmonic functions of time: 

Mj 

x(n) = Z [Aj exp(inoj) + Aj* exp(-inOj*)] , n = 1,2,. . .,N. (7) 

j=i 

That is, the AR time series (1) (with a(n)=0) and the series (7) are equivalent through equation (6). 
Equation (6), here referred to as the Prony’s relation, plays the central role in the present study, and 
is proved in the Appendix. It has been proposed as the basis of a method for harmonic analysis by 
Chao & Gilbert (1980). If we let o =u j+iOj, it is obvious from (7) that cOj is simply the angular fre- 
quency and ^ the exponential decay constant of the jth complex harmonic component in the time 
series. For this reason, Oj’s are called the complex frequencies. Their geometrical interpretation, 
according to Prony’s relation (6), is given in Figure 1 . The complex quantities j Aj ; j=l ,2,. . .,Mj j in 
(7) give the amplitude and phase of each complex harmonic component, and are called the complex 
amplitudes. They contain the same information as that in the M initial values of the AR time series 
(1) in the sense that the former can be determined by the latter, and vise versa. We shall see later 
that the last statement has important implications. 

Now let us study the case where some poles lie on the real axis of the z-plane. For example, if 
M is odd, we are bound to have at least one real pole. It is easy to show that real poles also satisfy 
Prony’s relation (in fact, the so-called Prony’s method was originally designed for fitting real expo- 
nentials, see Appendix); and they can be considered here as special cases of complex poles. Thus, a 
negative pole corresponds to a component which is an alternating real exponential series (a case of 
little practical interest), while a positive pole corresponds to a real exponential function of time (it 
could be a decaying exponential, a constant, or a growing exponential depending on whether the 
pole is less than, equal to, or greater than I, respectively). 
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Z-PLANE 



Figure 1 . Geometrical interpretation of the position of one complex conjugate pair of poles on the 
z-plane (denoted by®) in relation to the corresponding complex frequency a=cj+ioi. 


4 


Thus, for ME/AR spectral analysis the Prony’s relation in essence provides a connection between 
the “z-domain” representation and the time-domain representation of a random process. In the z 
domain a time series is modeled by Mj complex conjugate pairs of poles, whereas in the time domain 
it is modeled as a linear combination of Mj complex harmonic functions. (It follows that the AR 
linear prediction of a time series is in fact the extrapolation of these complex harmonic functions 
into the future.) Each complex conjugate pair of poles corresponds to a complex harmonic function, 
in a way summarized in Figure 1. In general, the number is much smaller than the number of data 
points N. By contrast, the conventional Fourier spectral analysis decomposes a time series into N/2 
real harmonic components (or pure sinusoids) at the fixed frequency interval 2ir/N. [In this sense, 
the Fourier analysis “models” a time series (signal as well as noise) exactly, in the form of the inverse 
Fourier transform.] This “brute-force” quantization of frequencies is, of course, what is responsible 
for the low resolving power of the Fourier analysis. The “data-adaptive” (Lacoss, 1971) ME/AR 
scheme does not suffer from this limitation and essentially gives a high-resolution “line” spectrum, 
as we shall see in the next section. 

The Meaning of the ME/AR Spectrum 

The function P(z) (equation 4) is a positive-valued function defined in the complex z-plane (or 
the z domain); the value of P(z) at the poles grows to infinity as an inverse-square function. The 
ME/AR spectrum P(o>) (equation 5) is the “cross section” of P(z) along the unit circle Cq (or the 
frequency domain, which is a one-dimensional subset of the z domain). It can be easily shown from 
equation (4) or (5) that P(oo) is inversely proportional to the square of the distance between the 
point exp(ico) (on Cq) and any of the poles. Therefore any pole that is close to Cg will give rise to a 
conspicuous peak in the spectrum— the closer the pole to C 0 the higher the peak, as depicted sche- 
matically in Figure 2. 

Now let us study the spectrum in the vicinity of a pole which lies close to C 0 and is isolated on 
the z-plane in the sense that the position and shape of the associated peak are not affected appreciably 
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by the presence of other poles. Let the pole be exp(ia 0 ) where ffQ=cj 0 +iG! 0 (see Figure 2); and we 
know from above that this pole (along with its complex conjugate) corresponds to a complex harmon- 
ic component with frequency coq and decay constant cv 0 . The fact that exp(ioQ) lies close to Cq im- 
plies that |a 0 |«l, or that the complex harmonic function has a high quality factor, Qs=co 0 /2a 0 »l. 
Thus, we have P(ou)=c<a 2 >/d 2 , where d is the distance between exp(ico) and exp(ia 0 ), and c is essen- 
tially a constant of a> depending on the distances between exp(ico) and all other poles, which vary 
little for co in the vicinity of coq. Invoking the cosine law and Taylor expansion, we get the following 
expression for the ME/AR spectrum of a high-Q, complex harmonic function in the vicinity of its 
frequency: 


P(co) = 


c<a 2 > 

“o + ( w " w 0 ) 2 


( 8 ) 


Equation (8) is a bell-shaped function of co with maximum at co=coq, as shown in Figure 3. This is 
essentially the basis for the usage of ME/AR spectrum in harmonic analysis of time series. 


Simple as it is, equation (8) has some important implications: 

(i) In general, the numerator <a 2 > does contain information about the amplitude of the com- 
plex harmonic function. Despite this, however, the spectrum P(ce) is largely determined by 
the denominator, or, the pole configuration— a conspicuous peak merely indicates the pre- 
sence of a pole lying close to Cq. Therefore in practice there is little information about the 
signal amplitude in the ME/AR spectrum. This is indeed a logical derivative of an earlier 
statement that the information on the complex amplitude is contained in the initial values 
of the AR model, not in the AR parameters. Therefore it is not surprising that the ME/AR 
spectral analysis under some circumstances may introduce spurious spectral peaks that does 
not correspond to any “signal”— they are simply the consequence of spurious poles (due to 
noises even though the noise level may be low) which happen to be close to Cq. 

(ii) The shape of the spectral peak is characterized by A, its full spectral width at half power 
(see Figure 3). From equation (8) it is immediately seen that A=2|a 0 i. Hence it is possible 
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to estimate from A the decay constant a 0 > anc * hence Q, of the complex harmonic compo- 
nent, in precisely the same way as would be done in the Fourier spectral analysis. This pro- 
cedure has been employed by various investigators in studying the Earth’s polar motion 
(Currie, 1974, Vincente & Currie, 1976, Graber, 1976) although, as far as I know, no theo- 
retical justification were previously given. However, since equation (8) depends on the 
square of a 0 , it cannot distinguish between positive and negative a 0 . In other words, a de- 
caying harmonic function (with a pole inside Cq) and a growing harmonic function (with a 
pole outside Cq) will have identical spectrum as long as they have the same |a 0 |, as first 
pointed out by Chao (1983). 

Lacoss, in a pioneering work (Lacoss, 1971), claimed that the ME/AR spectral height is propor- 
tional to the square of the signal-to-noise ratio (S/N). This is inconsistent with equation (8) and (i) 
above. Lacoss (1971) also claimed that the spectral width is inversely proportional to S/N, so that 
the integrated area under a spectral peak is proportional to S/N. This, again, is inconsistent with our 
argument (ii) above. In fact, according to our equation (8) the spectral area is 2<a 2 >/|a 0 |, which, 
incidentally, grows unbounded as the decay constant a 0 approaches zero (as for a pure sinusoid). At 
any rate, the term “power density spectrum” seems inappropriate for ME/AR spectrum because it 
does not reflect the nature of the latter. 

One (perhaps provocative) point to be made here is the following. Since the spectrum P(co) is 
only the “cross section” of the function P(z) in the frequency domain, it alone does not tell the whole 
story. Indeed, P(oo) is nothing but one convenient, but ambiguous way to display visually the content 
of a time series. What is important is the pole configuration of equation (2). In practice all the poles, 
or the zeros of the polynominal H(z) (equation 3), can be found numerically by means of, say, Bair- 
stow’s method (see e.g. Wilkinson, 1965). Bairstow’s method looks for real quadratic factors of H(z) 
one by one, which can then be further factorized into a complex conjugate pair of poles as the right 
side of equation (7). The latter, according to Prony’s relation, corresponds to a complex harmonic 
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function whose complex frequency can thus be determined. For example, a pole lying close to Cq 
(with absolute value «1) can be associated with a quasi-sinus oi dal component without looking at the 
ME/AR spectrum P(co). Moreover, it becomes pointless to discuss the interference of two spectral 
peaks in the frequency domain as long as their poles can be determined. 

Conclusions 

(1) The poles of a real autoregressive (AR) process are either real or in complex conjugate 
pairs. 

(2) Each complex conjugate pair of poles in the z domain corresponds to one complex harmon- 
ic function in the time domain. Special cases include pure sinusoids and real exponentials; 
the latter correspond to real poles. The relation between the position of the poles and the 
corresponding complex frequencies is provided by the Prony’s relation. Therefore, an AR 
model is one that models in a least-squares sense the time series as a linear combination -of 
complex harmonic function of time. 

(3) The complete information of the AR parameters of a random process is contained in the 
pole configuration in the (two-dimensional) z domain. The maximum-entropy /autoregres- 
sive (ME/AR) spectrum, being defined on the (one-dimensional) frequency domain, is only 
a convenient, but ambiguous visual representation. 

(4) Tile position and shape of an ME/AR spectral peak is determined by the corresponding 
complex frequency (which, in turn, is determined by the pole position relative to the Unit 
circle), and the height or area of the spectral peak contains little information about the 
amplitude of the harmonic components. The latter is contained, instead, in the initial 
values of the time series. 
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APPENDIX 


PRONY’S RELATION 

The so-called Prony’s method, due to Prony (1795), was originally designed for analyzing a 
linear combination of (real) exponential functions (see e.g. Froberg, 1969). It can be easily extended 
to encompass complex harmonic functions, as we shall do presently. The notations used are the same 
as in the text. 

Suppose we have the following equally spaced, real-valued time series: 

Mj 

x(n) = 2 [Aj exp(in<jj) + Aj* exp(-in<7j*)] , n = 1 ,2,. . .,N. (Al) 

Let exp (iaj)=Cj. Then 

x(l) = 2 [Aj Cj + Aj* Cj*] 
j 

x(2) = 2 [Aj Cj 2 + Aj* Cj 2 ] (A2) 

j 

x(M) = 2 [Aj Cj M + Aj * Cj * M ] 
j 

where M=2Mj . Now form the polynomial of degree M 

Mj 

H(z) = n (z — Cj) (z — Cj*) (A3) 

j=l 

= z M -Sj z M -1 - . . ,-S M . 

Note that the coefficients 1, -Sj ,. . .,— Sjyj are real. Then for any n>M, 

x(n) - Sj x (n-1) - S 2 x (n-2) - ... - Sj^ x (n-M) 

Mj 

= 2 Aj Cj n-M [cM- Sj CjM- 1 - . . . - Sjfl] + [compl. conj.] 

j=l 

Mj 

= 2 Aj Cj n_M H(Cj) + [compl. conj.] 

j=l 

?=0 [because H(Cj)=0] 

A-l 


(A4) 



or, 


M 

x(n) = 2 S k x (n-k), n = M+l . .,N. (A5) 

k=l 

Thus, equation (A3) relates the linear combination (Al) to the AR time series (A5). It is called the 
Prony’s relation. 
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